How Many Iterations in the Gibbs Sampler ?

When the Gibbs sampler is used to estimate posterior distributions (Gelfand and Smith, 1990), the question of how many iterations are required is central to its implementation. When interest focuses on quantiles of functionals of the posterior distribution , we describe an easily-implemented method for determining the total number of iterations required, and also the number of initial iterations that should be discarded to allow for \burn-in". The method uses only the Gibbs iterates themselves, and does not, for example, require external speciication of characteristics of the posterior density. Here the method is described for the situation where one long run is generated, but it can also be easily applied if there are several runs from diierent starting points. It also applies more generally to Markov chain Monte Carlo schemes other than the Gibbs sampler. It can also be used when several quantiles are to be estimated, when the quantities of interest are probabilities rather than full posterior distributions, and when the draws from the posterior distribution are required to be approximately independent. The method is applied to several diierent posterior distributions. These include a multivariate normal posterior distribution with independent parameters, a bimodal distribution, a \cigar-shaped" multivariate normal distribution in ten dimensions, and a highly complex 190-dimensional posterior distribution arising in spatial statistics. In each case the method appears to give satisfactory results. The results suggest that reasonable accuracy may often be achieved with 5,000 iterations or less; this can frequently be reduced to less than 1,000 if the posterior tails are known to be light. However, there are frequent \exceptions" when the required number of iterations is much higher. One important such exception is when there are high posterior correlations between the parameters; even crude correlation-removing The authors are grateful to Jeremy York for providing the data for Examples 4 and 5, for helping with the analysis and for useful discussions and suggestions, and to Julian Besag and an anonymous referee for helpful comments. A Fortran program called \Gibbsit" that implements the methods described here may be obtained from StatLib by sending an e-mail message to [email protected] containing the single line \send gibbsit from general". While the program is not maintained, questions about it may be addressed by e-mail to Adrian Raftery at [email protected]. 1 reparameterizations can greatly increase eeciency in such cases. Another important exception arises in hierarchical models when the Gibbs sampler tends to get \stuck"; …